function [F,X,X1,X2,X3,X4,Y1,Y2,Y3,Y4] = piecewise_polynomial_approximation(input_width,lut_width)
    step = 2^-input_width; 
    
    % data generation
    X = 1:step:2-step;
    F = cos(pi/2.*log2(X));
    F_deriv = (-pi*sin((pi*log(X))/(log(4))))./(X*log(4));
    
    % segment ranges
    seg_len = length(X)/4;
    seg1 = seg_len*0+1:seg_len*1;
    seg2 = seg_len*1+1:seg_len*2;
    seg3 = seg_len*2+1:seg_len*3;
    seg4 = seg_len*3+1:seg_len*4;
    
    % segment processing
    X1 = X(seg1);
    X2 = X(seg2);
    X3 = X(seg3);
    X4 = X(seg4);
    F1 = F(seg1);
    F2 = F(seg2);
    F3 = F(seg3);
    F4 = F(seg4);
    F1_deriv = F_deriv(seg1);
    F2_deriv = F_deriv(seg2);
    F3_deriv = F_deriv(seg3);
    F4_deriv = F_deriv(seg4);
    Y1 = linear_approx(X1,F1,F1_deriv,lut_width);
    Y2 = linear_approx(X2,F2,F2_deriv,lut_width);
    Y3 = linear_approx(X3,F3,F3_deriv,lut_width);
    Y4 = linear_approx(X4,F4,F4_deriv,lut_width);
end

function [Y1 s p] = linear_approx(X1,A1,A1_deriv,lut_width)
    len = length(X1);

    a_start = A1(1);
    a_end = A1(len);
    x_start = X1(1);
    x_end = X1(len);
    s = (a_end-a_start)/(x_end-x_start);
    s = floor(s*2^lut_width)*2^-lut_width;
    
    p_range = 1:0.00001:3;
    p_range_len = length(p_range);
    E = zeros(1,p_range_len);
    for i = 1:p_range_len,
        p = p_range(i);
        Y1 = p+s*(X1-a_start);
        e = find_max_error(Y1,A1);
        E(i) = e;
    end
    
    e_min = min(E);
    p_e_min_i = find(E == e_min);
    p = p_range(p_e_min_i);
    p = floor(p*2^lut_width)*2^-lut_width;
    Y1 = p+s*(X1-a_start);
end

function [e,e_i] = find_max_error(Y1,A1)
    e = 0;
    e_i = 0;
    len = length(Y1);
    for i = 1:len,
        e_t = abs(A1(i)-Y1(i));
        if(e_t > e)
            e = e_t;
            e_i = i;
        end
    end
end

